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Abstract 

Attitude estimation is often more difficult for spinning spacecraft than for three-axis stabilized platforms due to the 
need to follow rapidly- varying state vector elements and the lack of three-axis rate measurements from gyros. The 
estimation problem simplifies when torques are negligible and nutation has damped out, but the general case 
requires a sequential filter with dynamics propagation. This paper describes the implementation and test results for 
an extended Kalman filter for spinning spacecraft attitude and rate estimation based on a novel set of variables 
suggested in a paper by Markley [AAS 93-330] (referred to hereafter as Markley variables). 

Markley has demonstrated that the new set of variables provides a superior parameterization for numerical integ- 
ration of the attitude dynamics for spinning or momentum-biased spacecraft. The advantage is that the Markley var- 
iables have fewer rapidly- varying elements than other representations such as the attitude quaternion and rate vector. 
A filter based on these variables was expected to show improved performance due to the more accurate numerical 
state propagation. However, for a variety of test cases, it has been found that the new filter, as currently imple- 
mented, does not perform significantly better than a quaternion-based filter that was developed and tested in parallel. 

This paper reviews the mathematical background for a filter based on Markley variables. It also describes some 
features of the implementation and presents test results. The test cases are based on a mission using magnetometer 
and Sun sensor data and gyro measurements on two axes normal to the spin axis. The orbit and attitude scenarios 
and spacecraft parameters are modeled after one of the THEMIS (Time History of Events and Macroscale Interac- 
tions during Substorms) probes. 

Several tests are presented that demonstrate the filter accuracy and convergence properties. The tests include torque- 
free motion with various nutation angles, large constant-torque attitude slews, sensor misalignments, large initial 
attitude and rate errors, and cases with low data frequency. It is found that the convergence is rapid, the radius of 
convergence is large, and the results are reasonably accurate even in the presence of unmodeled perturbations. 

INTRODUCTION 

There are many ways to represent spacecraft attitude and its time rate of change. One representation has been 
put forward as superior in some respects for spinning or momentum-biased spacecraft. Markley (Ref. 1) has pre- 
sented a new set of variables that have fewer rapidly varying elements for spinning spacecraft than other commonly 
used representations. For this reason they provide numerical advantages when integrating the equations of motion. 
An earlier paper (Ref. 2) presented the mathematical background for an extended Kalman filter (EKF) designed to 
estimate spacecraft attitude using these new variables. The current paper draws heavily from Ref. 2 for a review of 
that background, with corrections, and then describes new work on the filter implementation and test results. 

Attitude estimation methods for non-spinning, three-axis stabilized spacecraft often make use of the Euler 
symmetric parameters (commonly called the quaternion) to represent the attitude (Ref. 3) with respect to an inertial 
reference frame. The quaternion is a four-component, globally nonsingular attitude representation. The state vector 
typically has seven components consisting of the quaternion plus three elements that provide bias corrections to the 
measured rates. Spinning spacecraft usually do not carry three-axis rate-sensing gyros, so the rotation rate vector 
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itself must be estimated rather than just the biases. Rates are needed for numerical time -propagation since the 
differential equation describing the attitude dynamics is 2 nd -order in time. A disadvantage of this state representation 
is that all four quaternion components and two of the body frame rate components are varying even for the simple 
case of torque-free nutating motion. As shown in Ref. 1 , alternative descriptions can be devised that have as few as 
three varying elements for this case, and thus, are easier to integrate numerically. These alternative representations 
will be referred to here as Markley variables. 

The goal of the current work is to design and build a Kalman filter for attitude estimation using the Markley 
variables. This new filter has been fully developed and incorporated into the attitude ground support system used for 
spinning spacecraft at the NASA Goddard Space Flight Center. In addition, a second new filter has been developed 
for comparison. This second filter is a more conventional EKF based on the quaternion and body rotation rate rather 
than Markley variables. It is found that the two filters have very similar performance; the Markley variable filter, as 
currently implemented, is not significantly more accurate than the quaternion and rate filter. 

The Markley variables are defined in Section 2. Section 3 reviews the partial derivatives required for a Kalman 
filter based on these variables. In particular, the observation partials are needed for the sensitivity matrix, and the 
linearized dynamics matrix is needed for time-propagation of the covariance matrix. Section 4 describes how the 
constraint condition can be explicitly incorporated in the filter equations, yielding a reduction of the filter update and 
the covariance matrix from 7 to 6 dimensions. Section 5 presents several numerical tests using simulated data 
patterned after the THEMIS spacecraft. Comparisons are made with a more typical EKF based on the quaternion 
and body rotation rate. Section 6 gives some conclusions. 

MARKLEY VARIABLES 

The seven variables of the attitude dynamics representation described in Ref. 1 are: 

• the angular momentum, L h in an inertial reference frame (typically, the geocentric inertial (GCI) frame), 

• the angular momentum, L s , in the body frame, 

• and a rotation angle, £ defined below. 

(Bold characters will be used throughout to denote 3-vectors.) A constraint condition holds for these variables: the 
magnitude of the angular momentum vector must be the same in the inertial and instantaneous body frames. 

The angle £ can be defined in a number of different ways. In this work, Markley’ s second representation is 
used, as follows. An intermediate frame is defined by rotating the inertial frame angular momentum into the body 
frame angular momentum, treating these as distinct vectors rather than as two representations of the same abstract 
vector. (Or, saying this more carefully, one transforms Lj into a frame where its components happen to be the same 
as those of L a .) The usual attitude matrix is only one of an infinite set of possible transformations that carry Lj into 
L b . The intermediate frame here is defined as the one obtained by choosing the most direct path from Lj to L B . For 
this purpose, the rotation vector is in the direction of the cross product L/ x L B . Call this transformation from the 
inertial frame to the intermediate frame A EI . Then, it can be shown that the transformation from the intermediate 
frame to the body, A BEi is simply a rotation about L B , as follows. 

The attitude matrix, A Bh is the transformation from the inertial to the body frame; thus, 



(i) 

By construction of A&, it is also true that 


Ajr/L/ — L b * 

(2) 

The transformation, A BEy is defined such that A BE A EI = A B i> Hence, one can write 


Afi£A £/ Ly = A fl/ Ly , 

(3) 

and this becomes 


A be L b — L b . 

(4) 


i 
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Thus, L b is an eigenvector of the orthogonal matrix A be- Equivalently, L B is the Euler rotation axis for the transfor- 
mation A be . The corresponding Euler rotation angle defines the 7 th Markley variable, f With this definition, A BE can 
be written (Ref. 3) 

. y r 1 cosC T t t sinC r 1 

A BE =c os^I-\ J 2 L b L b j—[L b x\, ( 5 ) 

where / is the 3x3 identity matrix, and the cross-product matrix is defined for any vector v to be 

0 -v, ■ (6) 

v, 0 _ 

(Properties of the cross-product matrix are extensively described in Ref. 3.) The angle ^is the rotation about L B to 
get from the intermediate frame to the body frame, and thus, represents predominantly the spin phase. 


The equation of motion for 


( 7 ) 


( 8 ) 

where 

(o BI =J-'{L B -h B ) (9) 

is the body rotation rate vector, J is the spacecraft moment of inertia tensor, and h B is the angular momentum 
relative to the body frame of any internal moving components. The external torque expressed in the body frame is 

N B ~ A B iN [ . 

Reference 1 shows that the time -dependence of the phase angle, £ satisfies the equation 

dC L [{L b +L i )-oj b ,+ L~ 2 (L b xL,)-(N b +N,)] 
dt L 2 + L B ■ L, 

KALMAN FILTER 

This section discusses the main parts of the Kalman filter. These are the propagation step and the update step for 
both the state vector and the state covariance matrix. Other details about Kalman filtering can be found, for example, 
in Ref. 4, with specific application to spacecraft attitude quaternion estimation in Ref. 5. 

Time Propagation 

The equations of motion can be integrated numerically to propagate the full 7-component state forward in time 
between sensor observations. As currently implemented, the integrator uses a 4th-order Runge-Kutta method, with 



Equations of Motion 

The Kalman filter requires expressions for the time-dependence of the state vector, 
the angular momentum in the inertial frame is 

dt 1 

where N t is the total external torque. This equation expressed in the body frame is 

~r~ = - a> B/ X L B , 

dt 


0 - 
M= 

r v > 


3 


r 


options to express the total torque as a sum of the command torques and environmental torques from gravity 
gradient and residual magnetization perturbations. 

The state error covariance is defined as the expectation value of the outer product of the errors in the state, 
P = E[{X - X ) (X - X ) r |, where X is the expectation value of the state. The covariance P is a 7x7 matrix; however, 
since the norms of L B and Lj are constrained to be equal, P is of rank 6 for an optimal estimator. Anticipating this 
problem and the resolution that will be given below, the full state vector now is defined to be 


X s 


<r 


( 12 ) 


where L ; is the unit vector in the direction of the inertial frame angular momentum. This state still has seven com- 
ponents, but the constraint now is that the first three elements have unit norm. 

The 7x7 error covariance satisfies the equation 

r!P 

— = FP + PF T +Q(t) , (13) 

dt 


where Q(t) is the process noise (Ref 4), and F is obtained from the state dynamics equation, 


with 


dt 


Fs V(xj) 


ax 


(14) 


(15) 


evaluated at the current state estimate. The function flX,t) is given by Eqs. 7, 8, and 11, except Eq. 7 now becomes 




dt 


(16) 


to account for the unitized inertial frame angular momentum in Eq. 12. With this change, the expression for F is 



(l, ■ N,I + l,nJ )/l 

-(i-l,l t ,)n,l t b/l 3 

^3x1 

F = 

®3x3 

[L B x}r'-[j-'L B x] 

^3x1 


d,L/(L 2 + L,-L b ) 

d B LI{l} +L, ■ L B ) 

Oixl 


(17) 


with 


and 


d, =-L T B ^f + l{j-'l b J +(n b + Nj [ l b x] 

dt 


d B = -(4 + L,J^- + (. r'Lj - (N s + N , ) r [L, x]/L 

dt 


(18) 


+(L b + L,) t J-' +(L, •y'Lgfc 
-[L b xL,) {N b +N i )L t b /L ■ 


( 19 ) 
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Sensitivity Matrix 


The sensitivity, H , is defined as the matrix of partials of the sensor observation with respect to the state X. This 
matrix is constructed in this section by chaining together partial derivatives with respect to intermediate variables. 
The current filter implementation is designed to accept observations from vector sensors, quaternion sensors (such 
as autonomous star trackers), or rate measurements from 1-, 2-, or 3-axis gyros. This section discusses H for vector 
and rate measurements. 

Vector observations, for example, may be obtained from a three-axis magnetometer or constructed from the 
output of an Earth sensor or a V-slit Sun or star sensor. It is assumed the vector observations have been transformed 
to the body frame where they are modeled as 

v? = A bi v, +n 

- exp( [a x] )A al v, + n (20) 

= (/ + [«x])v B 

= V B -[ V B x]« , 


where v 7 and v B are the reference vector expressed in the inertial frame and body frame, respectively, A est is the 
current estimate of the attitude matrix, a is the negative of the unknown attitude error vector, and n represents 
random noise with covariance R . Thus, the partial of v B obs with respect to a is -[v fi x]. Using a in the error state 
rather than -or simplifies signs elsewhere when constructing a quaternion-based filter. 

The partial of a with respect to the quaternion is related to the S matrix in Ref. 5 and can be shown to be 


3a 

a* 


<?4 

<73 

"?2 

-?i 


“9j 

<74 

9i 

-<7 2 

(21) 

. ^2 

-<7i 

<74 

“9j. 



Next, the attitude quaternion can be expressed as a function of the Markley variables, X, and the partials of the 
vector part of the quaternion, q, with respect to X can be determined. This yields 


dq_ 

at 


C, cos(£72) , C 2 sin(£/2) 


Lyj2(L 2 +L r L B ) V2(t 2 + t,-tJ ’ 


( 22 ) 


dq _ C 3 cos(£ 72 ) C 4 sin(^/ 2 ) 

dL B i(l 2 71, • l b ) ^i{l 2 + l, ■ l b ) ’ 


and 


where 


a? _ (t fi xt ; )sin(^/2) | (Lg + L ; )cos(£ 72 ) > 
a£ 2L^2(l 2 + L, ■ L b ) 2^2 (l 2 + L, L b ) ’ 


c 3 


, _ (L b xL,)LL t J 2 r 1 

L 2 + L~L~ {b] 


, (L 8 +L ; )LL r B /2 

L 2 + L r L B 


(L b xL,^L b+ ^ + ^lL b 
l } + Lf ‘ L B 


-[L,x], 


(24) 


(25) 

(26) 


(27) 
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and 


C 4 =/ + 


L,L T B (Lfl + L; ( Lfl+ 2 + 2L 2 ' Lb 


L 2 + L r L b 


(28) 


The partial of the ^-component of q is 


&h 1 

M 2 V2 {L 2 + L, ■ Lj 


COSl 


{C/2)L t b , cos (<r / 2)# (/ - L b L t b ) , - sin(£/2) — 


(29) 


The 3x7 sensitivity matrix for vector observations is obtained by combining these pieces using the chain rule. 


H = 


ax 


( da ) 

dq 

dq 

dq 


_d£, ’ 

dL B ' 



(30) 


This expression is evaluated approximately using the current state estimate on the right-hand side. 

Note that Eqs. 19 through 30 have been corrected and differ slightly from Ref. 2. All the partial derivatives have 
now been validated by extensive numerical tests. 

For rate observations, the actual measurements are the projections of the true spacecraft rotation rate along the 
sensitive axis of each gyro unit. (It is assumed that any bias, scale factor, or misalignment corrections have already 
been applied. Calibration methods are a topic of ongoing work.) If there are three sensitive axes, aligned with the 
unit vectors U u U 2 , and £/ 3 , then the sensor frame rate observations are simply I/j-fifci plus noise, for /=! to 3. As 
indicated by Eq. 9, the partial derivative with respect to the Markley variables is nonzero only for L s , thus 



Vx 

U 2 

U 3 


'^Bl 
' ft *8/ 


— IPjx3 U J ^3x1 . ’ 


(31) 


where 


Uj = 


C 

Ul 

\ul 


(32) 


If there are fewer than three gyro axes, H is obtained by deleting the appropriate rows from Eq. 3 1 . 


REDUCED REPRESENTATION 

Knowledge of the 7x7 state error covariance derives from an a priori estimate plus updates based on sensor 
observations. This information should account for the constraint condition at each step. If this is the case, the 
covariance matrix will have rank 6. Any treatment where the covariance has rank 7 must be using the initialization 
or update information sub-optimally. 

It can be difficult numerically to maintain the covariance P as a rank 6 matrix during the filter propagation and 
update steps. To avoid this problem entirely, it is preferable to recast the filter in terms of a 6-component error state, 
reducing the dimensionality of the covariance matrix and the filter update. 

The state reduction depends on two facts. First, as given in Eq. 12, the first three components of the state are 
taken to be the unit vector in the direction of L/ rather than Lj itself. Second, the filter update step will affect only the 
components perpendicular to the a priori angular momentum; the parallel component is affected only in 2 nd -order 
because it is a unit vector. Thus, any reference to the component parallel to L/ can be discarded, reducing the 
dimensionality of the state by one. 
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The rest of this section gives the details needed to define the reduced state and to use it to create a reduced 
representation of the filter. 


Unit Vector Covariance 


The 3x3 submatrix of P corresponding to the uncertainty in the unit vector angular momentum can be related to 
the covariance corresponding to the 3-vector L } as follows. Let the angular momentum and its error be written 


L ' = L + 6L , 

where SL has zero mean. Then, the corresponding unit vector and its error are 

L'={L+ SLfL 1 +26L-L+SL 1 )' ,/2 
.L + (l-LL')^ , 

where the subscript / has been dropped here for simplicity. The angular momentum covariance is 

p L = e[slsl t ] , 

and, from Eq. 34, the covariance of the unitized angular momentum vector is 

p l m l- 2 (i-LL t )p l (i-LL t ) ■ 


(33) 


(34) 


(35) 


(36) 


This property will be used below. 


Reduction to 6-Component State 


Next, construct the matrix M r that rotates Lj to be parallel to the Z-axis of the inertial frame. With this 
transformation, it is possible to reduce the dimensionality of the state simply by discarding the third component of 
M T L t , which is manifestly equal to unity. Only the X- and Y- components are kept since they can deviate from zero 
in the 1 st order at the update step of the filter. 

To this end, define 

M = 7cosz?+(l~cosz?)jw r +sin^[wx] , (37) 


where is the angle between the Z-axis and L h and 


ZxL 

n =7-7 — 

ZxL, 

With this definition, one has 

MZ = L l , 


by construction. Also, define the matrices 






1 

o' 

'l 

0 

o' 

, m 32 = 

0 

1 

0 

1 

0 







0 

0 


that actually perform the change of dimensionality. One can show that 

j A A™, 

M m, 2 m 2J M =I-L,L, ■ 


(38) 


(39) 


(40) 


(41) 
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Now, define the 6-component state to be 


x = 


m 23 M T L } 


= S T X 


(42) 


where S is the 6x7 matrix 


S T = 


m 23 M 0 2x3 0^ 


I 0 


®lx3 ^1x3 


(43) 


The m 2 3 matrix in Eq. 42 explicitly removes the Z-component of the rotated unit angular momentum vector. 
Equation 42 also shows that 


3x 

ax 


= . 


(44) 


Two important properties of the S matrix are that S T S - / 6x6 , and, using Eq. 41, 


i^ = ss t +z l . 


(45) 


In these expressions, / 6x 6 and / 7x7 are the 6x6 and 7x7 identity matrices, and the 7x7 matrix Z L is defined as 


L,L] 

®3x3 

®3xl 

^3x3 

®3x3 

®3xl 

O 

S 

®lx3 

®lxl 


(46) 


Reduced Form for the Filter 

The 6x6 covariance matrix for the reduced state now can be defined as 

P=S T PS (47) 

(overbar indicates 6-component reduced form). From Eq. 13, the time derivative is 

dt dt 

= s t f{ss t +z l )ps ( 48 ) 

+ S T p(sS r +zjF r S + S r G(OS , 

where the identity / 7x7 has been introduced in the form given in Eq. 45. The time-dependence of S has been 
neglected. In the absence of external torques, S is exactly constant, and in practice, the variation is small over each 
short numerical integration interval. Recomputing S once for each time step should be sufficiently accurate for 
propagation of the covariance. 

Note that 


PZ L =Z L P = 0 


(49) 
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since L f in Z L is annihilated by the unit vector covariance, as given by Eq. 36. Thus, if one defines 


f = s t fs 


( 50 ) 


and 

q=s t qs , 


(51) 


then Eq. 48 reduces to 


— = FP + PF r 
dt 


+Q(t) ■ 


(52) 


This is of the same standard form as Eq. 13 and carries the same information content but is not subject to any 
constraint condition. 

Similarly, the 7x3 Kalman gain 

K = PH T (HPH T + r) ] (53) 

can be written in a reduced 6x3 form by introducing / 7x7 expressed as in Eq. 45. Thus, 

s t k = s t p{ss t + z l )h t [h{ss t +z l )p{ss t + z l )h t + r\' 

= PH t [hPH t + r\' (54) 

= K , 


where R is the covariance of the sensor noise, n, in Eq. 20, and where the reduced 3x6 sensitivity matrix is 


H s HS 


(55) 


Finally, the state update is 

Sx = K(v?-v?). (56) 

The 3 rd , 4 th , and 5 th elements of the update vector Sx are added to the a priori estimate of L B . The 6 th element of 3x is 
added to the a priori estimate of £ The 1 st and 2 nd elements of Sx are corrections to the unitized angular momentum 
vector in the rotated frame; the corresponding a priori elements of x are zero by construction. Thus, the updated 
inertial frame angular momentum unit vector is 


Lj = M 


Sx { 

Sx, 


(l + Sx* + Sx \ ) 


(57) 


where M is given by Eq. 37. The full vector L I can be recovered at any time from L f and the norm of L B . 

In summary, the product SS T equals the identity plus a term that is annihilated by the projection operators in 
Eq. 36. Because of this property, the dimensionality of all the filter equations can be reduced using S and S T without 
changing the basic formalism. That is, the Kalman filter equations for the unconstrained 6-component state are 
formally identical to those for the original 7-component state. Only one extra step (Eq. 57) is required at each update 
to convert the 6-component error state into the full 7-component state that is needed for the next propagation step. 
The state error covariance is updated and propagated always in its 6x6 reduced form. 
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IMPLEMENTATION 


The spinning spacecraft EKF has been implemented in MATLAB as a subsystem of the Multimission Spin- 
Axis Stabilized Spacecraft (MSASS) Attitude Ground Support System (AGSS). The MSASS AGSS has been used 
for spinning spacecraft mission support for many years at the NASA/Goddard Space Flight Center (GSFC). The new 
EKF subsystem adds an important capability to solve for a time-dependent attitude history and could be used for 
real-time applications, if needed. 

The software is divided into an MSASS driver, an EKF main routine, and subroutines for time-propagation, 
computation of the sensitivity matrices, and filter update. The driver collects sensor measurement data and presents 
it to the filter as vector observations (rate and quaternion sensor data are not yet implemented in MSASS but can be 
generated in a simulator developed specifically to test the filter). If the commanded control torques are available, 
these should also be passed to the filter either as magnetic moments or directly as magnetic or thruster torques. 

After discarding any points flagged as bad, the EKF main routine checks the observation time steps and will 
interpose extra propagation steps, as needed, to improve the time-integration between sensor observations. The main 
routine also obtains the spacecraft ephemeris and geomagnetic field at each integration time step in order to compute 
environmental torques, if this option in the dynamics model is enabled. 

The environmental model includes torques due to gravity gradients and due to any residual constant magneti- 
zation of the spacecraft (see, e.g., Ref. 6). The gravity gradient torque can be written 

N GC =^fx(jf), (58) 

r 

where r is the spacecraft position vector, r is its magnitude, r is the corresponding unit vector, J is the spacecraft 
body frame inertia tensor, and ju is the Earth’s gravitational constant. The magnetization torque is simply 

=lxl(r 7 (mxB) , (59) 

where torque is in Newton-meters (N-m), m is any residual spacecraft magnetic moment in Ampere-meter 2 (A-m 2 ), 
and B is the Earth’s magnetic induction in milliGauss (mG), 

For each observation, the main routine calls the propagation subroutine to integrate the state vector and its 
covariance to the observation time. This subroutine uses a 4 th -order Runge-Kutta time-integration method. Next, the 
sensor residual and the appropriate sensitivity matrix are computed, and the state and covariance are updated. 

The initial state vector and covariance matrix are input to the EKF using familiar attitude and rate variables. The 
conversion to Markley variables is performed internally. Similarly, the estimated state and covariance are converted 
to attitude and rate variables before output. The internal use of Markley variables is totally transparent to the User. 

Singularity Avoidance 

It is clear throughout the development of this filter that the spacecraft angular momentum must be nonzero, but 
even with this stipulation, there are two other singularities. First, it can be seen from Eq. 1 1 that the g variable is 
undefined when L b -Lj = -L 2 . That is, the L B and Lj vectors cannot be 180 deg apart. For example, the angular 
momentum in the body frame typically is near the Z-axis. If this is so, this singularity precludes any attitude where 
the spin axis is at the celestial south pole. The second singularity can be seen in Eq. 38, which is needed to define 
the reduced representation. Problems occur when the inertial frame angular momentum is near ±Z, that is, near 
either celestial pole, regardless of where the angular momentum is in the body frame. 

It is an important feature of the software that it checks and maintains itself sufficiently far from these singular 
points of the representation. Whenever the estimated attitude is closer to any singular point than a user-specified 
tolerance, the software redefines the inertial reference frame. In this modified inertial frame, L r is far from the 
celestial poles and from -L B . The software transforms all reference vectors along with L h The inertial frame may be 
modified repeatedly whenever the solution approaches a singular point. Finally, to reconstruct the usual attitude, it is 
only a matter of bookkeeping to keep track of the times these extra rotations were applied. These transformations are 
all handled internally and again are totally transparent to the User. 
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TEST RESULTS 


Simulation Parameters 


A variety of tests have been performed using simulated data to exercise all the key features of the filter. The 
simulation is based on parameters from the THEMIS series of spinning spacecraft. For this mission, the simulation 
epoch is taken to be January 2, 2006 at 00 hours, 40 minutes, and the orbital parameters at epoch are 

semi-major axis = 43052.445 (km) 

eccentricity = 0.777800 

inclination = 0.16580627893946 (rad) 

right asc. of asc. node = 6.28177128263078 (rad) 

argument of perigee = 0.00276973678203 (rad) 

mean anomaly = 6.27689555620237 (rad) 


The initial true anomaly is 355.4 deg. 

The initial spin axis is directed 10 deg anti-sunward from the south ecliptic pole; the initial attitude quaternion is 
taken to be 


q Q = [-0.74356587177477; 

0.60417508447320; 

0.19420186426192; 

0.21063688554197] . (60) 


The spin rate is 20 revolutions per minute (rpm) approximately about the Z-axis. Various nutation angle offsets from 
the body Z-axis have been assumed for the tests. The moment of inertia tensor (kg-m 2 ) for a fully deployed 
configuration for THEMIS (Ref. 7) has been simplified to be axially symmetric. 


= 


200 

0 

0 


0 0 
200 0 
0 384 


(61) 


This simplifies numerical tests where the truth model is derived from closed-form analytic solutions to the dynamics 
equations. Other tests without axial symmetry and with nonzero off-diagonal terms in the inertia tensor have not 
caused problems. 

The sensor complement for THEMIS is a three-axis magnetometer (TAM) running at 8 observations per second 
(24 per spin period) and a slit Sun sensor that detects a Sun pulse once per spin and the Sun elevation angle from the 
nominal spin plane. For the EKF, the Sun elevation measurement and known slit azimuth angle are combined into 
an observed Sun unit vector. Gyros measure rotation rates about the two axes transverse to the spin during some 
parts of the mission. 

Unless otherwise indicated, the sensor noise in both the truth model and for EKF tuning are 0.1 deg for the Sun 
sensor elevation, 0.0007 sec for the Sun pulse time, and 1 mG per axis for the magnetometer. These are la values. 
The magnetometer noise covariance is the square of the ratio of this 1 mG error to the magnitude of the ambient 
geomagnetic field. The field is approximately 90 mG for the shorter tests and ranges from 90 to 10 mG for the 
longest runs (note that the orbit is quite eccentric). A rate error of 10‘ 4 deg/sec 1/2 is used for both the truth model and 
the EKF gyro noise R matrix. 

The process noise is taken, somewhat arbitrarily, to be the diagonal matrix 

Q(t) = diag^O" 6 , 10" 4 , 3X10" 4 , 1(T 7 , 10' 7 , 3 x 1 0 ” 7 ) , (62) 

corresponding to an error growth per hour of roughly 3.4 deg transverse to the spin direction, 5 deg around the spin 
axis, 1.1 deg/sec transverse rate, and 1.9 deg/sec spin rate. Note that this process noise is transformed on input from 
attitude and rate noise into the reduced Markley variable representation. 
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Test 1 - Torque-Free Rotation with Nutation 


To justify the added complexity, it must be shown that the new filter performs better than a more conventional 
quaternion-based filter. There are similarities between the new filter and quaternion-based filters (Ref 5). In particu- 
lar, the attitude and rate (or gyro rate bias) form a 7-element state vector subject to one constraint condition. In a 
quaternion-based filter, the constraint is that the quaternion must be normalized. Several of the tests presented here 
compare results from the new filter and the “unit vector filter” (UVF) (Ref. 8). The UVF has been used for attitude 
ground support at NASA/GSFC for many three-axis stabilized missions over the past 12 years. For application to 
spinning spacecraft, the UVF has been modified to estimate the quaternion and rotation rate rather than the gyro 
biases, to use dynamics propagation rather than gyro propagation, and to include a 7-dependent term in the 
linearized dynamics matrix to be used in Eq. 13: 


F = 


-M 

^3x3 


~7 

/-'([/wxHcux]/) 


(63) 


The first test scenario is simple torque-free motion with various values for the nutation angle. The filters were 
given 60 seconds convergence time, and statistics were collected over the next 90 seconds (30 spin periods). Table 1 
shows the mean error in the direction of the estimated angular momentum vector and the standard deviation of the 
error in the rotation phase about the body Z-axis. The new filter is noisier than the UVF in the spin phase estimate. 
At larger nutation angles, the new filter is slightly better at determining the angular momentum direction (which is 
operationally the more important statistic for most missions). However, the differences are small and the errors for 
both filters easily meet the 1 deg mission attitude knowledge requirement expected for THEMIS. 


Table 1. Comparison of New Markley Variable Filter and Unit Vector Filter for Torque-Free Nutation. 


Nutation 
Angle (deg) 

Mean Error in Ang. Momentum Direction (deg) 

Std. Dev. of Error about Body Z-axis (deg) 

New Filter 

UVF 

New Filter 

UVF 

0 

0.09 

0.09 

0.14 

0.10 

0.5 

0.09 

0.09 

0.15 

0.10 

1.0 

0.09 

0.09 

0.15 

0.10 

2.5 

0.08 

0.09 

0.14 

0.08 

5.0 

0.06 

0.08 

0.15 

0.07 

7.5 

0.07 

0.09 

0.18 

0.07 

10.0 

0.08 

0.10 

0.17 

0.06 


Test 2 - Attitude Reorientation 

The new filter behaves very well during large attitude reorientation slews. This test simulates a 90 deg attitude 
slew over a 230 sec time span. THEMIS is expected to perform large slews such as this during the early part of the 
mission prior to boom deployment. 

The filter followed the changing direction of the angular momentum vector throughout the slew. Figure 1 shows 
the attitude error relative to the truth model. The error in the estimated direction of the angular momentum, averaged 
over several runs, is only 0.09 deg during the slew. This is indistinguishable from the errors after slew completion. 
The errors for these tests using the UVF are 2 times larger, averaging 0.18 deg. 

Test 3 - Sensor Misalignments 

The TAM on THEMIS will be mounted on a deployable boom. The deployed alignment is expected to be repro- 
ducible only to within ±1 deg. To test how this would affect the attitude, the TAM misalignment was set to +1, +1, 
and -1 deg on the X-, F-, and Z-axes. Figure 2 shows a sample of the resulting attitude errors. The mean angular 
momentum pointing error is only 0.86 deg. 
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Attitude Errors Relative to Truth Model 
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Fig. 1. Attitude errors relative to truth model for a 90 deg slew. The slew ends at 230 seconds. 
The filter follows the slew with a mean pointing error of only 0.09 deg. 


Attitude Errors Relative to Truth Model 





Fig. 2. Attitude errors relative to truth model given a 1 deg misalignment on each TAM axis. 
The true attitude is modeled by torque-free motion with 5 deg nutation. 
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When the slit Sun sensor is also shifted by 1 deg in azimuth and 1 deg in elevation, the mean attitude errors are 
found to range between 0.5 to 1.9 deg, depending on the signs of the misalignment angles. The wide range occurs 
because the systematic errors from the Sun sensor misalignment can add to or subtract from those of the TAM. The 
errors with the UVF were similar. 

Test 4 - Sensor Complement 

This test investigates how performance varies with the data frequency of the TAM and with the availability of 
the two-axis gyro measurements. The attitude scenario is a 90 deg reorientation as in Test 2. 

Table 2 presents test results showing mean angular momentum pointing errors averaged over several runs. The 
pointing errors are slightly smaller for the new filter, but the difference is barely significant. What is more signifi- 
cant is that the convergence time for the UVF became longer and less predictable as the data frequency decreased. 
The new filter’s convergence time did not degrade in this way. 

The table also shows that using gyro data does not improve the attitude estimation. This may be due to the 
amount of gyro noise or may indicate filter tuning problems. The gyros do help the filter quickly to converge to the 
correct rotation rate when initialized with a large rate error (see Test 5). 

Perhaps the most interesting result indicated in Table 2 is how well the filter performs even when the TAM data 
frequency is so low there are fewer than one observation per spin period. 

Table 2. Comparison of New Markley Variable Filter and Unit Vector Filter for 90 deg Slew with Various Sensors. 


Sensor Complement 

Mean Error in Angular Momentum Direction During 
90 deg Attitude Reorientation Slew (deg) 

New Filter 

UVF 

Sun sensor, 8 Hz TAM, 4 Hz x- and y-axis gyro 

0.09 

0.18 

Sun sensor, 8 Hz TAM, no gyro 

0.09 

0.10 

Sun sensor, 1 Hz TAM, 4 Hz x- and y-axis gyro 

0.13 

0.20 

Sun sensor, 1 Hz TAM, no gyro 

0.13 

0.16 

Sun sensor, 1/3 Hz TAM, 4 Hz x- and y-axis gyro 

0.21 

0.24 

Sun sensor, 1/3 Hz TAM, no gyro 

0.17 

0.16 

Sun sensor, 0.1 Hz TAM, 4 Hz x- and y-axis gyro 

0.21 

0.28 

Sun sensor, 0.1 Hz TAM, no gyro 

0.27 

0.34 


The THEMIS orbit is very eccentric and has a large semimajor axis. The geomagnetic field is relatively weak 
over most of the orbit and will be useful as an attitude measurement only near perigee. For much of the orbit, the 
reference field model errors will be too large compared to the field magnitudes for TAM observations to provide 
attitude information. Two tests were performed to examine the effect of not having input from the TAM. First, a 
reorientation slew test was performed using only Sun and X,T-gyro data. The filter was allowed to converge using 
TAM data, but the TAM was disabled during most of the slew. The new filter again performed better than the UVF. 
The error in the angular momentum direction after the slew averaged 0.47 deg for the new filter and 0.85 deg for the 
UVF. In the next test, the filter was allowed to converge using Sun and TAM data for 200 seconds. The motion was 
torque-free with various nutation angles from zero to 2 deg. Solutions were obtained using only Sun sensor data for 
the next several thousand seconds. The error growth could easily be kept under a few degrees by choosing the 
dynamics integration step size. However, it was most surprising that the UVF performed much better for these 
torque-free Sun-only tests. For a 4000 sec test with 1 deg nutation angle, the UVF showed no error growth compared 
to Sun and TAM solutions, while the new filter showed linear error growth. 
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Test 5 - Convergence Tests 


The final series of tests examined how reliably the filter converges when initialized far from the true state. Over 
50 runs were performed with a wide sampling of large initial attitude errors. In most of these cases, the initial atti- 
tude errors were between 90 and 180 deg. Both the new filter and the UVF converged to within 0.3 deg of the true 
attitude in less than 2 minutes (40 spin periods). Most cases took less than 30 seconds to converge, with only a small 
number taking the full 2 minutes. 

When both the attitude and rate were initialized with large errors, the filters still usually converged quickly. 
Both filters converged within 2 minutes when the initial rotation rate error was 15 percent on the spin axis and 200 
percent on the X- and T-axes with 90 to 180 deg initial attitude errors. 

If gyro data is not included, there are some extreme cases where the new filter did not converge within 4 
minutes. In these tests, the initial X- and T-rate errors were 300 percent of the true rate and the initial attitude error 
was over 90 deg. With gyro data, even these difficult cases converged quickly. 

CONCLUSIONS 

The new spinning spacecraft EKF has been shown to perform well under several test scenarios. Its initial con- 
vergence is very robust and the overall accuracy is good. However, the filter does not perform significantly better 
than a more conventional quaternion-based filter such as the UVF. Additional tests and work on the software will 
look into whether improvements can be made that bring the new filter up to its original expectations. Regardless, 
this work has succeeded in providing two powerful filters for attitude determination. This adds a significant new 
capability to the attitude ground support system used for spinning spacecraft missions at the NASA Goddard Space 
Flight Center. 

Ongoing work is investigating how to determine sensor biases, scale factors, misalignments, and time-tag 
errors, either along with the attitude state or in a parallel utility. Several types of biases are already estimated as part 
of the current batch-method spinning spacecraft attitude ground support system, where the spin vector is assumed to 
be constant for the entire batch. However, it may be useful for early mission support to be able to determine the 
biases before the rotation has fully damped down to simple principal axis spin. 
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